Quick Start Guide

Get started with practical examples using GDAL Scripts

Overview

This guide walks you through three practical examples demonstrating common tasks with GDAL Scripts. Each example uses real scripts from the repository with actual command signatures.

These examples assume you have completed the installation steps and have test planetary imagery available. You can download sample data from USGS Astrogeology Data.

Example 1: Coordinate Conversion

One of the most common tasks in planetary data processing is converting between different coordinate systems. The gdal2Coordinates directory provides four utilities for this purpose.

Converting Pixel Coordinates to Latitude/Longitude

The pixel2longlat.py script reads the coordinate system and geotransformation matrix from an image and reports latitude/longitude coordinates for a specified pixel.

1. Understand the Script

   Script location: gdal2Coordinates/pixel2longlat.py

   Purpose: Given pixel (sample) and line coordinates, report the latitude/longitude coordinates for the center of that pixel.

   Usage:

   python pixel2longlat.py sample line infile
   

   Parameters:

   • sample: Pixel column number (x-coordinate)
   • line: Pixel row number (y-coordinate)
   • infile: Input geospatial image file
2. Run the Example

   Find the coordinates of pixel (100, 200) in a lunar image:

   cd gdal2Coordinates
   python pixel2longlat.py 100 200 lunar_mosaic.tif
   

   Expected Output:

   pixel: 100              line: 200
   longitude: -15.234567   latitude: 23.456789
   longitude: 15d14'4.44"W latitude: 23d27'24.44"N
   

   The script outputs coordinates in both decimal degrees and DMS (degrees, minutes, seconds) format.

3. Try Related Conversions

   The gdal2Coordinates directory includes complementary scripts:

   Lon/Lat to Meters

   python longlat2meters.py -15.234567 23.456789 lunar_mosaic.tif
   

   These scripts work together to convert between any coordinate system based on the image’s projection.

Understanding the Output

The coordinate conversion accounts for:

• The image’s geotransformation matrix (position, rotation, scale)
• The coordinate reference system (CRS) defined in the image
• Pixel center vs. pixel corner conventions

The script automatically shifts to the center of the pixel by adding half the pixel size to the calculated coordinates.

Example 2: Terrain Analysis with Slope Calculation

The gdal_baseline_slope.py script calculates specialized slopes using various baseline lengths, a method developed specifically for planetary terrain analysis.

1. Understand Baseline Slopes

   Script location: gdal_baseline_slope/gdal_baseline_slope.py

   Purpose: Calculate slopes using variable baseline lengths measured in pixels. Different baselines reveal terrain characteristics at different scales.

   Key Concepts:

   • Baseline = 1, 2, or 5: Uses corner pixels of a moving window
   • No baseline specified: Uses Horn’s Method (3×3 window, equivalent to gdaldem slope)

   Usage:

   python gdal_baseline_slope.py [-baseline VALUE] [-ot Byte] [-crop] infile outfile.tif
   
2. Calculate Standard Slope (Horn's Method)

   Without specifying a baseline, the script uses Horn’s Method:

   cd gdal_baseline_slope
   python gdal_baseline_slope.py mars_dem.tif mars_slope.tif
   

   Output:

   Warning: Using Horn's method for slope calculation, send -baseline VALUE to set specialized calculation.
   band: 1 complete.
   

   This produces a 32-bit floating-point slope raster in degrees.

3. Calculate Baseline Slope

   For specialized terrain analysis with a 5-pixel baseline:

   python gdal_baseline_slope.py -baseline 5 mars_dem.tif mars_slope_5baseline.tif
   

   The baseline parameter changes which pixels are used in the calculation:

   • Baseline 2: Uses 3×3 window corners
   • Baseline 5: Uses 6×6 window corners

   Larger baselines smooth out small-scale roughness and highlight regional slopes.

4. Generate 8-bit Output

   For visualization or to reduce file size, create 8-bit output:

   python gdal_baseline_slope.py -baseline 5 -ot Byte mars_dem.tif mars_slope_8bit.tif
   

   This scales slope values (0-50+ degrees) to 1-255 using:

   DN = (Slope * 5) + 0.2
   

   The script generates both 32-bit and 8-bit versions:

   • 32bit_mars_slope_8bit.tif (full precision)
   • mars_slope_8bit.tif (scaled to byte)
5. Optional: Crop Edge Pixels

   The -crop flag removes edge pixels affected by the baseline window:

   python gdal_baseline_slope.py -baseline 5 -crop mars_dem.tif mars_slope_cropped.tif
   

   This trims 1-5 pixels from edges depending on baseline value, ensuring all output pixels have complete neighborhoods.

Performance Considerations

The current implementation loads the full image into memory and can be slow for large datasets. Plan accordingly:

• Use -ot Byte to reduce memory usage for output
• Work with tiled or regional subsets for very large DEMs
• Expect processing time to scale with image size

Dependencies

This script requires:

• NumPy: Array operations
• SciPy: Generic filtering functions

conda install numpy scipy

Example 3: Clipping Data to Valid Ranges

The gdal_clip2range.py script is essential for cleaning planetary data by setting out-of-range values to NoData.

1. Understand the Use Case

   Script location: gdal_clip2range/gdal_clip2range.py

   Purpose: Clip pixel values to a defined valid range. Values outside this range are set to NoData.

   Common applications:

   • Removing erroneous values from I/F (reflectance) data
   • Eliminating elevation outliers from DEMs
   • Cleaning up processed imagery

   Usage:

   python gdal_clip2range.py infile outfile.tif min_valid max_valid [output_noData]
   
2. Clip with Existing NoData Value

   Most planetary images have a NoData value defined. Preserve it:

   cd gdal_clip2range
   python gdal_clip2range.py lunar_dem.tif lunar_dem_clipped.tif -100.5 500.2
   

   What this does:

   • Sets pixels < -100.5 meters to NoData
   • Sets pixels > 500.2 meters to NoData
   • Maintains the original NoData value from the input file
   • Preserves all pixels within the valid range

   Output:

   input file: lunar_dem.tif
   min: -100.5
   max: 500.2
   NoData: -3.40282e+38
   file created: lunar_dem_clipped.tif
   
3. Set Valid Range for Reflectance Data

   For I/F (reflectance) images that should only contain values between 0 and 1:

   python gdal_clip2range.py mars_if.tif mars_if_clipped.tif 0 1 0
   

   Note the fifth parameter: Here we explicitly set the output NoData value to 0 because:

   • The input file may not have NoData defined
   • For reflectance data, 0 is a reasonable NoData value
   • All invalid pixels (< 0 or > 1) will be set to 0

   Only specify the optional output_noData parameter if the input file has no NoData value defined. Otherwise, omit it to preserve the original NoData value.

4. Verify the Results

   Check that clipping worked as expected:

   gdalinfo -stats mars_if_clipped.tif
   

   Look for:

   • Minimum and maximum values should fall within your specified range
   • NoData value should be set correctly
   • Statistics should reflect the clipped data

Technical Details

The script:

1. Opens the input file and creates a copy for output
2. Reads each band as a NumPy array
3. Applies the range mask: data[(data < minValue) | (data > maxValue)] = noData
4. Writes the modified array back to the output file
5. Clears existing statistics metadata (which would be invalid)

For more information and community discussion, see: OpenPlanetary: Crop Image DN Range

Next Steps

Explore More Tools

Now that you’ve tried these examples, explore other categories:

Format Conversion

Convert to ISIS3, PDS4, or generate point clouds:

• gdal2ISIS3/Astropedia_gdal2ISIS3.py
• PDS4gdal/isis3_to_pds4_LOLA_pvl.py
• gdal2PLY/gdal2PLY.py

Projection Tools

Work with IAU coordinate reference systems:

• OGC_IAU2000_WKT_v2/Source_Python/create_IAU2000_wkt_v3.py
• NewCenterLon_Equi/NewCenterLon_Equi.py

Metadata Operations

Extract and manipulate image metadata:

• gdal2metadata/gdal2metadata.py
• gdal_copylabel/gdal_copylabel.py

Spatial Data

Work with vectors and spatial features:

• ogr_footprintinit2shp/footprintinit2shp.py
• ogr_isisminer2shp/isisminer2shp.py

Batch Processing

For processing multiple files, create shell scripts or Python wrappers:

Example Batch Script

#!/bin/bash
# Process all TIF files in a directory

for dem in *.tif; do
    echo "Processing $dem..."
    python gdal_baseline_slope.py -baseline 5 -ot Byte "$dem" "slope_${dem}"
done

Getting Help

Each script directory contains a README.md with detailed usage instructions:

cat gdal_baseline_slope/README.md
cat gdal2Coordinates/README.md
cat OGC_IAU2000_WKT_v2/README.md

Run any script without arguments to see its usage information:

python gdal_clip2range.py

Common Workflows

Landing Site Analysis Workflow

1. Download DEM of the region of interest
2. Clip to valid range using gdal_clip2range.py
3. Calculate slopes at multiple baselines:
   python gdal_baseline_slope.py -baseline 1 dem.tif slope1.tif
   python gdal_baseline_slope.py -baseline 5 dem.tif slope5.tif
   
4. Extract coordinates of potential sites using pixel2longlat.py
5. Export to shapefile for analysis in GIS software

Image Mosaic Preparation

1. Match image extents using gdal_match_image_extents.py
2. Clip to valid ranges for consistent data quality
3. Convert projection if needed using IAU WKT definitions
4. Mosaic using standard GDAL tools like gdal_merge.py

Data Archive Conversion

1. Convert from ISIS3 cubes to GeoTIFF for web services
2. Generate PDS4 labels using PDS4gdal scripts for archiving
3. Create metadata extracts using gdal2metadata.py
4. Generate footprint shapefiles for inventory database

Report Issues or Contribute

Found a bug or have a feature request? Visit the GitHub repository to open an issue or contribute improvements.